notebook.community

Edit and run

Methanol-Synthese

$T_{ein}=493,15K$

p=50 bar

Bilanzen

Stoffbilanzen

$\dot n_i = \dot n_{i, 0} + \sum_{j}{\nu_ij \xi_j}$

Energiebilanz

$\begin{array}{ll} 0 &= \dot Q + (\dot n (\Delta H^\circ_{(T)}-\Delta H^\circ_0))_{ein}- (\dot n (\Delta H^\circ_{(T)}-\Delta H^\circ_0))_{aus} + \sum\limits_{j}{\xi_j (-\Delta Hr_j(T))}\\ &= 0 + (\dot n (\Delta H^\circ_{(T)}-\Delta H^\circ_0))_{ein}- (\dot n (\Delta H^\circ_{(T)}-\Delta H^\circ_0))_{aus} + \sum\limits_{j}{\xi_j (-\Delta Hr_j(T))}\\ \end{array}\\ $

Gleichgewichtskonstanten

$\begin{array}{ll} exp \left(- \frac{\Delta G_i}{R T} \right) &= K_p K_{\phi^{eq}} = K_x \prod\limits_{i} \left( \frac{p}{p^0}\right)^{\nu_i} K_{\phi^{eq}} \\ &=\prod\limits_{i} (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum\limits_{i} \nu_i}(n)^{-\sum\limits_{i} \nu_i} K_{\phi^{eq}}\end{array}$

$p^0 = 1 bar$

Idealer Gas, $K_{\phi^{eq}}=1$

Methode A) Geringe Veränderung der Reaktionsenthalpie mit der Temperatur

Van't Hoff, $\frac{d ln K}{dT} = -\frac{\Delta H}{R T^2} \sim \Rightarrow ln \left(\frac{K}{K'} \right) = -\frac{\Delta H^0}{R}\left(\frac{1}{T} - \frac{1}{T'} \right)$

$\begin{array}{ll} K_{(493,15K)} &= K_{(298,15K)} \times exp \left[-\frac{\Delta H^0}{R}\left(\frac{1}{493,15 K} - \frac{1}{298,15 K} \right)\right] \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$

Methode B) Wechselwirkung der Reaktionsenthalpie mit der Temperatur [SVNA]

$\Delta H^\circ = \Delta H_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT}$

$\Delta S^\circ = \Delta S_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$

$\Delta G^\circ = \Delta H^\circ - T \Delta S^\circ = \Delta H_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} - T \Delta S_0^\circ - R T \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$

$\Delta S_0^\circ = \frac{\Delta H_0^\circ - \Delta G_0^\circ}{T_0}$

$\Delta G^\circ = \Delta H_0^\circ - \frac{T}{T_0}(\Delta H_0^\circ -\Delta G_0^\circ) + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} - R T \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$

$\begin{array}{ll} K_{(T)} &= exp \left(-\frac{\Delta H_0^\circ}{R T} + \frac{(\Delta H_0^\circ -\Delta G_0^\circ)}{R T_0} - \frac{1}{T}\int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} + \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}\right) \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$

Somit läßt sich K(T) bestimmen, insofern man über einen Ausdruck für $Cp_i(T)$ verfügt. Bei geringer Veränderung der Wärmekapazität Cp im Temperatur-Bereich kann man auch einen bestimmten Mittelwert als ~konstant einsetzen.

Methode C) Gibbs'sche Energie-Funktion - Gef(T) - aus Thermochemischen Tabellen [BP]

$-Gef(T) = -[G(T)-H(298,15)]/T$

$-R ln(K) = \sum\nu_i Gef_i - \sum \nu_i H_i(298,15K)/T$

In thermochemischen Tabellen [BP] sind die Werte -Gef(T) verfügbar.

Literaturhinweise

  • [SVNA] Smith J.M., Van Ness H.C., Abbott M.M.; Introduction to chemical engineering thermodynamics; 6th ed.; McGraw-Hill; New York; 2001; S. 458-462.
  • [BP] Barin Isan, Platzki Gregor; Thermochemical data of pure substances; 3. ed.; VCH; New York; 1995.

Adiabatischer Fall





In [1]:



    


from scipy import optimize
import numpy as np

p = 50. # bar
temp = 273.15 + 220. # K
t0_ref = 298.15 # K
r = 8.314 # J/(mol K)

namen = ['CO', 'H2', 'CO2', 'H2O', 'CH3OH']

n0co = 50.
n0h2 = 100.
n0co2 = 0.
n0h2o = 0.
n0ch3oh = 0.

ne = np.array([n0co, n0h2, n0co2, n0h2o, n0ch3oh])

nuij = np.array([[-1, -2, 0, 0, +1] ,
                 [0, -3, -1, +1, +1], 
                 [-1, +1, +1, -1, 0]]).T

h_298 = np.array(
    [-110.541, 0., -393.505, -241.826,-201.167]) * 1000 # J/mol

g_298 = np.array(
    [-169.474, -38.962, -457.240, -298.164, -272.667]) * 1000 # J/mol

# Berechne delta Cp(T) mit Temperaturfunktionen für ideale Gase (SVN).

# Koeffizienten für Cp(T)/R = A + B*T + C*T^2 + D*T^-2, T[=]K
# Nach rechts hin: A, B, C, D
# Nach unten hin: CO, H2, CO2, H2O, CH3OH
cp_coefs =  np.array([
    [
        y.replace(',', '.') for y in x.split('\t')
    ] for x in """
3,3760E+00	5,5700E-04	0,0000E+00	-3,1000E+03
3,2490E+00	4,2200E-04	0,0000E+00	8,3000E+03
5,4570E+00	1,0450E-03	0,0000E+00	-1,1570E+05
3,4700E+00	1,4500E-03	0,0000E+00	1,2100E+04
2,2110E+00	1,2216E-02	-3,4500E-06	0,0000E+00
""".split('\n') if len(x)>0], dtype=float)

def cp(t):
    return r * (
        cp_coefs[:,0] + 
        cp_coefs[:,1] * t + 
        cp_coefs[:,2] * t**2 + 
        cp_coefs[:,3] * t**-2
    ) # J/(mol K)

# Berechne H(T), G(T) und K(T) mit Cp(T)

def h(t):
    return (
        h_298 + 
        r * cp_coefs[:,0]*(t-t0_ref) + 
        r * cp_coefs[:,1]/2.*(t**2-t0_ref**2) + 
        r * cp_coefs[:,2]/3.*(t**3-t0_ref**3) -
        r * cp_coefs[:,3]*(1/t-1/t0_ref)
    ) # J/mol

def g(t, h_t):
    return (
        h_t - t/t0_ref*(h_298 - g_298) -
        r * cp_coefs[:,0]*t*np.log(t/t0_ref) -
        r * cp_coefs[:,1]*t**2*(1-t0_ref/t) - 
        r * cp_coefs[:,2]/2.*t**3*(1-(t0_ref/t)**2) +
        r * cp_coefs[:,3]/2.*1/t*(1-(t/t0_ref)**2)
    ) # J/mol

def k(t, g_t):
    delta_g_t = nuij.T.dot(g_t)
    return np.exp(-delta_g_t/(r * t))


delta_gr_298 = nuij.T.dot(g_298)

delta_hr_298 = nuij.T.dot(h_298)

cp_493 = cp(493.15) # J/(mol K)
h_493 = h(493.15) # J/mol
g_493 = g(493.15, h_493) # J/mol
k_493 = k(493.15, g_493) # []

for i, f in enumerate(delta_hr_298):
    print('Delta H_' + str(i+1) + '(298.15K)=' + str(f/1000.) + 'kJ/mol')

print('\n')
for i, f in enumerate(k_493):
    print('K' + str(i+1) + '(493K)=' + str(f))
print('\n')

n0 = np.array([n0co, n0h2, n0ch3oh])
    
def fun(x_vec):    
    nco = x_vec[0]
    nh2 = x_vec[1]
    nch3oh = x_vec[2]
    xi1 = x_vec[3]
    t = x_vec[4]
    
    n = np.array([nco, nh2, nch3oh])
    
    cp_t = cp(t)
    h_t = h(t)
    g_t = g(t, h_t)
    k_t = k(t, g_t)
    
    h_ein = h_493[[1, 2, -1]]
    cp_ein = cp_493[[1, 2, -1]]
    cp_t = cp_t[[1, 2, -1]]
    h_t = h_t[[1, 2, -1]]
    g_t = g_t[[1, 2, -1]]
    
    delta_h_t = nuij[[1, 2, -1]].T.dot(h_t) # J/mol
    
    f1 = -nco + n0co - xi1
    f2 = -nh2  + n0h2 -2*xi1
    f3 = -nch3oh + n0ch3oh +xi1
    f4 = -k_t[0] * (nco * nh2**2) + \
            nch3oh * (p/1.)**-2 * (nco + nh2 + nch3oh)**-(-2)
    #f5 = np.sum(
    #    np.multiply(n0, h_ein) - 
    #    np.multiply(n, h_t)
    #) + xi1 * (-delta_h_t[0])
    f5 = np.sum(
        np.multiply(n0, cp_ein)*493.15 - 
        np.multiply(n, cp_t)*t 
    ) + xi1 * (-delta_h_t[0])
    
    return [f1, f2, f3, f4, f5]

x0 = np.append(n0, [0., 493.15])

sol = optimize.root(fun, x0)
f_final = - sol.x[:3].reshape([3,1]) + ne[[0,1,4]].reshape([3,1]) + nuij[:,0][[0,1,4]].reshape([3,1])*sol.x[-2]

print(sol)
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)

print('\n\n')
print('T_ein=493.15K, p=50 bar, in adiabatischem Reaktor')
print('Lösung für nur einzige Reaktion (ohne CO2):\n')
for i, f in enumerate(sol.x[:2]):
    print('n_' + namen[i] + '= ' + str(f) + ' mol')
print('n_' + namen[-1] + '= ' + str(sol.x[2]) + ' mol')
print('T= ' + str(sol.x[-1]) + ' K')

n0 = np.array([n0co, n0h2, n0co2, n0h2o, n0ch3oh])
n0 = ne

# Lösung des einfacheren Falls in schwierigerem Fall einwenden.
def fun(x_vec):    
    nco = x_vec[0]
    nh2 = x_vec[1]
    nco2 = x_vec[2]
    nh2o = x_vec[3]
    nch3oh = x_vec[4]
    xi1 = x_vec[5]
    xi2 = x_vec[6]
    xi3 = x_vec[7]
    t = x_vec[8]
    
    n = np.array([nco, nh2, nco2, nh2o, nch3oh])
    xi = np.array([xi1, xi2, xi3])
    
    h_ein = h_493
    cp_ein = cp_493
    cp_t = cp(t)
    h_t = h(t)
    g_t = g(t, h_t)
    k_t = k(t, g_t)
    
    delta_h_t = nuij.T.dot(h_t) # J/mol
    
    f1 = -nco + n0co - xi1 +0 -xi3
    f2 = -nh2  + n0h2 -2*xi1 -3*xi2 +xi3
    f3 = -nco2 + n0co2 +0 -xi2 +xi3
    f4 = -nh2o + n0h2o +0 +xi2 -xi3
    f5 = -nch3oh + n0ch3oh +xi1 +xi2 -0
    f6 = -k_t[0] * (nco * nh2**2) + \
        nch3oh * (p/1.)**-2 * (nco + nh2 + nco2 + nh2o + nch3oh)**-(-2)
    f7 = -k_t[1] * (nco2 * nh2**3) + \
        nch3oh * nh2o * (p/1.)**-2 * (nco + nh2 + nco2 + nh2o + nch3oh)**-(-2)
    f8 = -k_t[2] * (nco * nh2o) + \
        nco2 * nh2 * (p/1.)**0 * (nco + nh2 + nco2 + nh2o + nch3oh)**-0
    f9 = np.sum(
        np.multiply(n0, (h_ein-h_298)) - 
        np.multiply(n, (h_t-h_298))) + np.dot(xi, -delta_h_t)
    #f9 = np.sum(
    #    np.multiply(n0, cp_ein)*493.15 - 
    #    np.multiply(n, cp_t)*t) + np.dot(xi, -delta_h_t)
    
    return [f1, f2, f3, f4, f5, f6, f7, f8, f9]

x0 = np.append(n0, [0., 0., 0., sol.x[-1]])

sol = optimize.root(fun, x0)

f_final = - sol.x[:5].reshape([5,1]) + ne.reshape([5,1]) + nuij.dot(sol.x[5:-1].reshape([3,1]))

print('\n\n')
print('T_ein=493.15K, p=50 bar, in adiabatischem Reaktor.')
print('Lösung für alle drei Reaktionen, mit CO2:\n')
for i, f in enumerate(sol.x[:5]):
    print('n_' + namen[i] + '= ' + str(f) + ' mol')

print('\n')

for i, f in enumerate(sol.x[5:-1]):
    print('xi_' + str(i) + '= ' + str(f) + ' mol')
    
print('\n')
    
print('T=' + str(sol.x[-1]) + ' K, oder...')
print('T=' + str(sol.x[-1]-273.15) + ' °C')

print('\n')
print('0 = Q + Sum(Delta H)_ein - Sum(Delta H)_aus')
bilanz = np.sum(
    np.multiply(n0, (h_493-h_298)) -
    np.multiply(sol.x[:5], (h(sol.x[-1])-h_298))
) + np.dot(sol.x[5:-1], -nuij.T.dot(h(sol.x[-1]))) 
annaeherung = np.sum(
    np.multiply(n0, cp_493)*493.15 -
    np.multiply(sol.x[:5], cp(sol.x[-1]))*sol.x[-1]
) + np.dot(sol.x[5:-1], -nuij.T.dot(h(sol.x[-1])))
print('-Q = (n.(H_t-H_298))_ein -(n.(H_t-H_298))_aus  + Sum(xi_j * (-Delta Hr_j)) = ' + 
      str(bilanz) + 'J/h')
print('-Q = (n Cp T)_ein - (n Cp T)_aus + Sum(xi_j * (-Delta H_j)) = ' + 
      str(annaeherung) + 'J/h' + 
      '(nur Annäherung. Fehler: ' + '{:.2f}'.format(
          (annaeherung-bilanz)/bilanz) + ')' )
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)



    


















    






Delta H_1(298.15K)=-90.626kJ/mol
Delta H_2(298.15K)=-49.488kJ/mol
Delta H_3(298.15K)=-41.138kJ/mol


K1(493K)=0.0088102868389
K2(493K)=5.7133657404e-05
K3(493K)=154.204845956


    fjac: array([[ -6.98220403e-05,  -2.13470694e-16,  -9.52744512e-21,
         -5.12921178e-03,  -9.99986843e-01],
       [  1.57166023e-02,  -9.97133966e-03,  -4.60048510e-18,
          9.99813608e-01,  -5.12942058e-03],
       [ -3.73103145e-01,  -6.50699281e-01,   6.61350208e-01,
         -6.24396114e-04,   2.92538676e-05],
       [ -8.20521478e-01,   5.63874061e-01,   9.19102575e-02,
          1.85216484e-02,  -3.77114697e-05],
       [  4.32771916e-01,   5.08465399e-01,   7.44424883e-01,
         -1.73205518e-03,  -2.13332210e-05]])
     fun: array([ -4.44089210e-16,  -8.88178420e-16,   0.00000000e+00,
        -5.52020651e-10,   6.98491931e-10])
 message: 'The solution converged.'
    nfev: 20
     qtf: array([  1.44182179e-06,  -1.56560613e-06,   9.40392430e-10,
        -2.90851842e-08,   2.75554739e-09])
       r: array([  1.43221252e+04,   2.25742130e+04,   3.02905388e+04,
        -2.03411166e+05,   7.23437740e+03,   1.00287427e+02,
         9.86723160e+01,  -3.95568925e+03,   1.27484710e+02,
        -1.51205819e+00,   1.01056495e+01,  -2.68089360e-01,
        -6.18361166e+01,   1.94702164e+00,  -2.21688703e-03])
  status: 1
 success: True
       x: array([  46.22169757,   92.44339514,    3.77830243,    3.77830243,
        614.06573205])



Zustand der Optimisierungs-Funktionen

[[ -4.44089210e-16]
 [ -8.88178420e-16]
 [  0.00000000e+00]]



T_ein=493.15K, p=50 bar, in adiabatischem Reaktor
Lösung für nur einzige Reaktion (ohne CO2):

n_CO= 46.2216975716 mol
n_H2= 92.4433951433 mol
n_CH3OH= 3.77830242836 mol
T= 614.065732048 K



T_ein=493.15K, p=50 bar, in adiabatischem Reaktor.
Lösung für alle drei Reaktionen, mit CO2:

n_CO= 45.3678348003 mol
n_H2= 90.7356696005 mol
n_CO2= 1.85405940675e-16 mol
n_H2O= 1.4430304831e-17 mol
n_CH3OH= 4.63216519974 mol


xi_0= -2.64802822417 mol
xi_1= 7.28019342391 mol
xi_2= 7.28019342391 mol


T=606.759730614 K, oder...
T=333.609730614 °C


0 = Q + Sum(Delta H)_ein - Sum(Delta H)_aus
-Q = (n.(H_t-H_298))_ein -(n.(H_t-H_298))_aus  + Sum(xi_j * (-Delta Hr_j)) = -1.21246557683e-07J/h
-Q = (n Cp T)_ein - (n Cp T)_aus + Sum(xi_j * (-Delta H_j)) = -18389.5024134J/h(nur Annäherung. Fehler: 151670305242.10)



Zustand der Optimisierungs-Funktionen

[[  0.00000000e+00]
 [  0.00000000e+00]
 [ -1.85405941e-16]
 [ -1.44303048e-17]
 [  0.00000000e+00]]

















In [3]:



    


print('Lösung, in 30 Dezimalzahlen')
print('')
for part in sol.x:
    print('{:.30g}'.format(part).replace('.',','))



    


















    






Lösung, in 30 Dezimalzahlen

45,367834800264532191249600146
90,7356696005290643824992002919
1,85405940674631041373621166986e-16
1,44303048310322819158516687663e-17
4,63216519973546780875039985403
-2,64802822417277949895719757478
7,28019342390824775179680727888
7,28019342390824775179680727888
606,759730614496675116242840886

Min(G)

$\Delta G_{f i}^{0} + R T ln(y_i \hat{\phi_i} P/P^0) + \sum\limits_{k}{ \lambda_k a_{ik}}=0 \hspace{10mm} (i=1,2,...,N)$

$\sum\limits_{i}{n_i a_{ik}}=A_k \hspace{10mm} (k = 1,2,...,\omega)$

Content source: santiago-salas-v/walas

Similar notebooks:

notebook.community | gallery | about